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Abstract. Coleman's theory of p-adic integration figures prominently 
in several number-theoretic applications, such as finding torsion and ra- 
tional points on curves, and computing p-adic regulators in if-theory 
(including p-adic heights on elliptic curves). We describe an algorithm 
for computing Coleman integrals on hyperelliptic curves, and its imple- 
mentation in Sage. 



1 Introduction 

One of the fundamental difficulties of p-adic analysis is that the totally discon- 
nected topology of p-adic spaces makes it hard to introduce a meaningful form 
of antidifferentiation. It was originally discovered by Coleman that this problem 
can be circumvented using the principle of Frobenius equivariance. Using this 
idea, Coleman introduced a p-adic integration theory first on the projective line 
[9], then (partly jointly with de Shalit) on curves and abelian varieties [10], [8]. 
Alternative treatments have been given by Besser [3] using methods of p-adic 
cohomology, and by Berkovich [2] using the nonarchimedean Gel'fand transform. 

Although Coleman's construction is in principle quite suitable for machine 
computation, this had only been implemented previously in the genus case 
[5]. The purpose of this paper is to present an algorithm for computing single 
Coleman integrals on hyperelliptic curves of good reduction over C p for p > 2, 
based on the third author's algorithm for computing the Frobenius action on the 
de Rham cohomology of such curves [17]. We also describe an implementation 
of this algorithm in the Sage computer algebra system. 

For context, we indicate some of the many potential applications of explicit 
Coleman integration. Some of these will be treated, with additional numerical 
examples, in the first author's upcoming PhD thesis. (Some of these applications 
will require additional refinements of our implementation; see Section 5.) 



— Torsion points on curves. Coleman's original application of p-adic integration 
was to find torsion points on curves of genus greater than 1. This could 
potentially be made effective and automatic. 



— p-adic heights on curves. Investigations into p-adic analogues of the con- 
jecture of Birch and Swinnerton-Dyer for Jacobians of hyperelliptic curves 
require computation of the Coleman-Gross height pairing [11]. This global 
p-adic height pairing can, in turn, be decomposed into a sum of local height 
pairings at each prime. In particular, for C a hyperelliptic curve over Q p 
with p a prime of good reduction and for D\,D 2 6 Div°(C) with disjoint 
support, the Coleman-Gross p-adic height pairing at p is given in terms of 
the Coleman integral [10] 



for an appropriately constructed differential lud 1 associated to the divisor 
D\. This pairing is effectively computable by work of the first author [1]. 
Using this work, it should be possible (using ideas of Besser [4]) to add in 
local heights away fromp, and thus compute the Coleman-Gross height pair- 
ing on Jacobians of hyperelliptic curves. (In genus 1, one can then compare 
to an alternate computation based on work of Mazur-Stein-Tate [22] and 
Harvey [16].) 

— p-adic regulators. A related topic to the previous one is the computation 
of p-adic regulators in higher X-theory of arithmetic schemes, which are 
expected to relate to special values of L-functions. Some computations in 
genus have been made by Besser and de Jeu [5]. 

— Rational points on curves: Chabauty 's method. For C a smooth proper curve 
over Z[-^], the Chabauty condition on C is that rank J(C) (Z [-^j ) < dim J(C), 
where J(C) denotes the Jacobian of the curve. When the Chabauty condi- 
tion holds, there exists a 1-form uj on J(C) an with J Q w = for all points 
P G J{C) (Z We might be able to compute C(Z[^]) if we can find all 

points P G C an such that J„ P u> — 0. This method has already been used in 
many cases, by Coleman and many others; see [23] for a survey (circa 2007). 
To apply Chabauty's method in a typical case, one needs the integral of to at 
some point in a residue disc, with which one can find all zeroes of the integral 
in the residue disc. Several methods are suggested in [23, Remark 8.3] for 
doing this, including Coleman integration. However, no serious attempt has 
been made to use numerical Coleman integration in Chabauty's method; it 
seems likely that it can handle cases where the other methods suggested in 
[23, Remark 8.3] for finding constants of integration prove to be impractical. 

— Rational points on curves: nonabelian Chabauty. It may be possible to use 
(iterated) Coleman integration to find rational points on curves failing the 
Chabauty condition, using Kim's nonabelian Chabauty method [18]. As a 
demonstration of the method, Kim [19] gives an explicit double integral 
which vanishes on the integral points of the minimal regular model of a 
genus 1 curve over Q of Mordell-Weil rank 1. The erratum to [19] includes a 
corrected formula, together with some numerical examples computed using 
the methods of this paper. 

— p-adic polylogarithms and multiple zeta values. These have been introduced 
recently by Furusho [13], but little numerical data exists so far. 
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2 Coleman's theory of p-adic integration 

In this section, we recall Coleman's p-adic integration theory (for single integrals 
only) in the case of curves with good reduction. This theory involves some con- 
cepts from rigid analytic geometry which it would be hopeless to introduce in 
such limited space; some standard references are [6] and [12]. (See also [10, §1].) 

Let C p be a completed algebraic closure of Q p , and let O be the valuation 
subring of C p . Choose once and for all a branch of the p-adic logarithm, i.e., a 
homomorphism Log : C* — > C p whose restriction to the disc {x G C p : \x — 1| < 
1} is given by the logarithm series log(x) = X^i(l~ x )V*- (The choice of branch 
has no effect on the integrals on differentials of the second kind, i.e., everywhere 
meromorphic differentials with all residues zero.) 

We first introduce integrals on discs and annuli within P 1 . 

Definition 1. Let I be an open subinterval of [0,+oo). Let A{I) denote the 
annulus (or disc) {t G : \t\ G /}. For J2iez dt G ^\^yc p and P,Q G 



This is easily shown not to depend on the choice of the coordinate t. 

Remark 2. Note that because of the division by i + 1 in the formula for the 
integral, we are unable to integrate on closed discs or annuli. 

We next turn to curves of good reduction. 

Definition 3. By a curve over O, we will mean a smooth proper connected 
scheme X over O of relative dimension 1. Equip the function field K(X) with 
the p-adic absolute value, so that the elements of K{X) of norm at most 1 
constitute the local ring in X of the generic point of the special fibre X of X. 

Let Xq denote the generic fibre of X as a rigid analytic space. There is a 
natural specialization map from Xq to X ; the inverse image of any point of X is 
a subspace of Xq isomorphic to an open unit disc. We call such a disc a residue 
disc of X. 



A(I), define 




Definition 4. Let X be a curve over O. By a wide open subspace of Xq, we 
will mean a rigid analytic subspace of Xq of the form {x <E Xq : \ f(x)\ < A} for 
some f € K{X) of absolute value 1 and some A > 1. 

Coleman made the surprising discovery that there is a well-behaved integra- 
tion theory on wide open subspaces of curves over O, exhibiting no phenomena 
of path dependence. (Note that one needs to consider wide open subspaces even 
to integrate differentials which are holomorphic or meromorphic on the entire 
curve.) In the case of hyperelliptic curves, Coleman's construction of these inte- 
grals using Frobenius lifts will be reflected in our technique for computing the 
integrals. For the general case, see [10, §2], [3, §4], or [2, Theorem 1.6.1]. 

Theorem 5 (Coleman). We may assign to each curve X over O and each 
wide open subspace W of Xq a map /lyy : Div (W) x f2^y/c — > C p , subject to 
the following conditions. (Here Div(T4 7 ) denotes the free group on the elements 
of W , and Div°(M / ) denotes the kernel of the degree map deg : Div(W) — > Z 
taking each element ofW to 1.) 

(a) (Linearity) The map fi\y is linear on Div°(W) and C p -linear on ^^y/ C • 

(b) (Compatibility) For any residue disc D of X and any isomorphism ip : W R 
D — > A(I) for some interval L, the restriction of nw 

to Div°(WnD)xf2l v/Cp 

is compatible with Definition 1 via if). 

(c) (Change of variables) Let X' be another curve overO, letW be a wide open 
subspace of X' , and let ip : W — >• W' be any morphism of rigid spaces relative 
to an automorphism ofC p . Then 

Hw>{1>(-),-)=liw{;P{-)). (1) 

(d) (Fundamental theorem of calculus) For any Q = J ^2 i Ci(Pi) S Div°(VF) and 
any f eO(W), » w (Q,df) ='E i c i f(P i ). 

Remark 6. One cannot expect path independence in the case of bad reduction. 
For instance, an elliptic curve over C p with bad reduction admits a Tate uni- 
formization, so its logarithm map has nonzero periods in general. In Bcrkovich's 
theory of integration, this occurs because the nonarchimedean analytic space 
associated to this curve X has nontrivial first homology. 



3 Explicit integrals for hyperelliptic curves 

We now specialize to the situation where p > 2 and X is a genus g hyperelliptic 
curve over an unramified extension K of Q p having good reduction. We will 
assume in addition that we have been given a model of X of the form y 2 = f{x) 
such that deg f{x) =2^ + 1 and / has no repeated roots modulo p. (This restric- 
tion is inherited from [17], where it is used to simplify the reduction procedure. 
One could reduce to this case after possibly replacing K by a larger unramified 
extension of Q p , by performing a linear fractional transformation in x to put 



one root at infinity, thus reducing the degree from 2g + 2 to 2g + 1 .) We will 
distinguish between Weierstrass and non-Weierstrass residue discs of X, which 
respectively correspond to Weierstrass and non-Weierstrass points of X. 

To discuss the differentials we will be integrating, we review a core definition 
from [17]. Let X 1 be the affine curve obtained by deleting the Weierstrass points 
from X, and let A = K[x, y, z]/ (y 2 — f(x), yz — 1) be the coordinate ring of X'. 

Definition 7. The Monsky- Washnitzer (MW) weak completion of A is the ring 
A' consisting of infinite sums of the form 

B i (x)eK[x],degB i <2g\, 




further subject to the condition that v p (Bi(x)) grows faster than a linear function 
of i as i — > ±oo. We make a ring out of these using the relation y 2 = f(x). 

These functions are holomorphic on wide opens, so we will integrate 1-forms 

w = s(z,2/)0, g(x,y)eAl (2) 

Note that we only consider 1-forms which are odd, i.e., which are negated by 
the hyperelliptic involution. Even 1-forms can be written in terms of x alone, 
and so can be integrated directly as in Definition 1. (This last statement would 
fail if we had taken A* to be the full p-adic completion of A, rather than the 
weak completion. This observation is the basis for Monsky- Washnitzer's formal 
cohomology, which is used in [17].) 

Note that the class of allowed forms includes those meromorphic differentials 
on X whose poles all belong to Weierstrass residue discs. For some applications 
(e.g., p-adic canonical heights), it is necessary to integrate meromorphic differ- 
entials with poles in non-Weierstrass residue discs. These will be discussed in 
[!]■ 

Note also that for ease of exposition, we describe all of our algorithms as if it 
were possible to compute exactly in A^ . This is not possible for two reasons: the 
elements of A^ correspond to infinite series, and the coefficients of these series 
are polynomials with p-adic coefficients. In practice, each computation will be 
made with suitable p-adic approximations of the truly desired quantities, so one 
must keep track of how much p-adic precision is needed in these estimates in 
order for the answers to bear a certain level of p-adic accuracy. We postpone 
this discussion to § 4.1. 



3.1 A basis for de Rham cohomology 

We first note that any odd differential co as in (2) can be written uniquely as 

co = df + c oj H lc 2g -\U2 g -\ (3) 




That is, the form a basis of the odd part of the de Rham cohomology of . 
The process of putting uj in the form (3), using the relations 



y 2 = /(*), 

d(zV) = (2^-y +1 +jx i f(x)yi- 1 ) ^, 



can be made algorithmic; see [17, §3]. (Briefly, one uses the first relation to 
reduce high powers of x, and the second to reduce large positive and negative 
powers of y.) Using properties from Theorem 5 (linearity and the fundamental 
theorem of calculus) , the integration of co reduces effectively to the integration 
of the u)i. 

It may be convenient for some purposes to use a different basis of de Rham 
cohomology. For instance, the basis x % dx/2y 3 (i = 0, . . . , 2g — 1) is crystalline 
(see the erratum to [17]), so Frobenius will act via a matrix with p-adically 
integral entries. 

3.2 Tiny integrals 

We refer to any Coleman integral of the form Jp u> in which P, Q lie in the same 
residue disc (Weierstrass or not) as a tiny integral. As an easy first case, we give 
an algorithm to compute tiny integrals of basis differentials. 

Algorithm 8 (Tiny Coleman integrals). 

Input: Points P,Q E X(C P ) in the same residue disc (neither equal to the point 
at infinity) and a basis differential Wj. 
Output: The integral Jp u>i. 

1. Construct a linear interpolation from P toQ. For instance, in a non- Weierstrass 
residue disc, we may take 



Remark 9. One can similarly integrate any u> holomorphic in the residue disc 
containing P and Q. If to is only meromorphic in the disc, but has no pole at 
P or Q, we can first make a polar decomposition, i.e., wrie w as a holomorphic 
differential on the disc plus some terms of the form cj(t — r) 1 , and integrate 
the latter terms directly. (If lo is everywhere meromorphic, this is achieved by a 
partial fractions decomposition.) 



x(t) = (1 - t)x(P) + tx(Q) 



y(t) = y/JWfi), 

where y{t) is expanded as a formal power series in t. 
2. Formally integrate the power series in t: 

rQ rQ a~ r 1 -vltY rlv( 




3.3 Non-Weierstrass discs 



We next compute integrals of the form J p u>i in which P,Q e ^(C p ) lie in dis- 
tinct non-Weierstrass residue discs. The method of tiny integrals is not available; 
we instead employ Dwork's principle of analytic continuation along Frobenius, 
in the form of Kedlaya's algorithm [17] for calculating the action of Frobenius 

on de Rham cohomology. Note that we calculate the integrals JpOJi for all i 
simultaneously (We modify the presentation in [17] by keeping track of exact 
differentials, which are irrelevant for computing zeta functions.) 

Algorithm 10 (Kedlaya's algorithm). 
Input: The basis differentials {wj} 2 ^ 1 . 

Output: Functions fa e and a 2g x 2g matrix M over K such that 4>*{wi) = 
dfi + J2j^o MijUJj for a p-power lift of Frobenius (p. 

1. Since K is an unramified extension ofQ p , it carries a unique automorphism 
4>k lifting the Frobenius automorphism x >-> x p on its residue field. Extend 
4>k to a Frobenius lift on A^ by setting 

<f>(x) = x p , 



i=0 



l/2\ (<j>K(f)(x p ) - fjxYf 
i J y 1 ^ 



noting the series converges in because <j>K(f)(x p ) — f(x) p has positive 
valuation. (This choice of cf){y) ensures that 4>{y) 2 — 4>(f(x)), so that the 
action on A^ is well-defined. 
2. Use a Newton iteration to compute y/4>(y). Then for i = 0, . . . , 2g— 1, proceed 
as in § 3.1 to write 

d 29-1 

0>o = p^- 1 ^^ = dfi + j2 M ^ (5) 



for some fa G and some 2g x 2g matrix M over K. 



We may use Algorithm 10 to compute Coleman integrals between endpoints 
in non-Weierstrass residue discs, as follows. (Note that our recipe is essentially 
Coleman's construction of the integrals in this case.) 

Algorithm 11 (Coleman integration in non-Weierstrass discs). 
Input: The basis differentials {wi} 2 ^ 1 , points P,Q E X(C P ) in non-Weierstrass 
residue discs, and a positive integer m such that the residue fields of P, Q are 
contained in F p ™ . 

Output: The integrals {JpUj}^ 1 . 



1. Calculate the action of the m-th power of Frobenius on each basis element 
(see Remark 12): 

2g-l 

(<j> m )*L0i = df i +J2 M ij"j- ( 6 ) 

3=0 

2. By change of variables (see Remark 13), we obtain 

29-1 ,Q r<t> m {P) rQ 

J2 (M - / Uj = MP) - MQ) - Ui- Ui (7) 

j=0 Jp Jp J<t> m {Q) 

(the fundamental linear system/ As the eigenvalues of the matrix M are 
algebraic integers of C p -norm p m / 2 ^ 1 (see [17, §2]), the matrix M — I is 
invertible, and we may solve (7) to obtain the integrals Jp Wj. 

Remark 12. To compute the action of </> m , first perform Algorithm 10 to write 

29-1 

(j)*LOi = dgi + ^ BijUj- 

3=0 

If we view /, g as column vectors and M, B as matrices, we then have 

/ = ^-Hg) + B<f> m - 2 (g) + ■■■ + B(f) K (B) ■ ■ ■ <t>T\B)g 
M = B4> K {B)---4>™-\B). 

Remark 13. We obtain (7) as follows. By change of variables, 

r<t> m (Q) 



r<P W) rW 

/ Ui= (<f> m )*"i 

J<f> m (P) JP 

rQ 29- 1 

= / (dfi + M ^ 
Jp 



3=0 

29-1 r Q 



MQ)-MP)+Y J Ma f u>j. 

Jp 



3=0 -> P 

Adding jp u>i + J^ m ^ to both sides of this equation yields 

rQ r<t> m (P) rQ 2 S -1 q 

">i= "i+ "i + MQ) - MP) + Y] Mij 

Jp Jp J<t> m (Q) j=0 Jp 

which is equivalent to (7). 
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Definition 14. A Teichmiiller point of Xq is a point fixed by some power of 
<p. Each non-Weierstrass residue disc contains a unique such point: if (x,y) € 



X is a non-Weierstrass point, the Teichmuller point in its residue disc has x- 
coordinate equal to the usual Teichmuller lift of x. This leaves two choices for 
the y-coordinate, exactly one of which has the correct reduction modulo p. Note 
that Teichmuller points are always defined over finite unramified extensions of 
Q P . 

Remark 15. A variant of Algorithm 11 is to first find the Teichmuller points 
P', Q' in the residue discs of P, Q, then note that from the fundamental linear 
system (7), we have 

2fl-l r Q> 

£(M-/) y / Uj = fi(P') - fi(Q'). (8) 
i=o Jp ' 

From (8), we obtain the integrals Jp, Wj. Finally, write Jp w, — Jp, uJi as the 
sum Jp uii + jS, LUi of tiny integrals. 



3.4 Weierstrass endpoints of integration 

Suppose now that P, Q lie in different residue discs, at least one of which is Weier- 
strass. Since a differential u> of the form (2) is not meromorphic over Weierstrass 
residue discs, we cannot always even define Jp to, let alone compute it. We will 
thus assume (to cover most cases arising in applications) that u> is everywhere 
meromorphic, with no pole at either P or Q. We then make the following obser- 
vation. 

Lemma 16. Letuj be an odd, everywhere meromorphic differential on X . Choose 
P,Q € X(C P ) which are not poles of u, with P Weierstrass. Then for i the hy- 

perelliptic involution, Jp uj = \J®q^u. In particular, if Q is also a Weierstrass 
point, then Jp uj = 0. 

Proof. Let I := Jp uj = Jp ®\—uj) = J^^uj. Then by additivity in the end- 
points, we have J®q^ w = 27, from which the result follows. 

If P belongs to a Weierstrass residue disc while Q does not, we find the 
Weierstrass point P' in the disc of P, then apply Lemma 16 to write 

,Q f P' 1 ,Q 

u= u + - u. (9) 

JP JP 1 Jl(Q) 

The first integral on the right side of (9) is tiny, while the second integral involves 
two points in non-Weierstrass residue discs, and so may be computed as in the 
previous section. The situation is even better if P, Q both belong to residue discs 
containing respective Weierstrass points P', Q': in this case, by Lemma 16, Jp uj 
equals the sum J p w + J uj of tiny integrals. 



Remark 17. Beware that Lemma 16 does not generalize to iterated integrals. For 
instance, for double integrals, if both integrands are odd, the total integrand is 
even, so the argument of Lemma 16 tells us nothing. It is thus worth considering 
alternate approaches for dealing with Weierstrass discs, which may generalize 
better to the iterated case. We concentrate on the case where P lies in a Weier- 
strass residue disc but Q does not, as we may reduce to this case by splitting 
j p uj = Jp ijj + j R uj for some auxiliary point R in a non- Weierstrass residue 
disc. 

In Algorithm 11, the form fi belongs to A> and so need not converge at P. 
However, it does converge at any point R near the boundary of the disc, i.e., in 
the complement of a certain smaller disc which can be bounded explicitly. We 
may thus write Jp uji — jp u>i + jp u>i for suitable R in the disc of P, to obtain 
an analogue of the fundamental linear system (7). Similarly, when we write u> as 
in (3), we can find R close enough to the boundary of the disc of P so that / 
converges at R, use (3) to evaluate j R uj, then compute jp uj as a tiny integral. 
One defect of this approach is that forcing R to be close to the boundary of the 
residue disc of P forces R to be defined over a highly ramified extension of Q p , 
over which computations are more expensive. 

An alternate approach exploits the fact that for P in the infinite residue 
disc but distinct from the point at infinity, we may compute jp uj directly using 
Algorithm 11. This works because both the Frobenius lift and the reduction 
process respect the subring of consisting of functions which are meromorphic 
at infinity. When P lies in a finite Weierstrass residue disc, we may reduce to 
the previous case using a change of variables on the rr-line to move P to the 
infinite disc. However, one still must use the approach of the previous paragraph 
to reduce evaluation of jp uj to evaluation of the jp uji. 

4 Implementation notes and precision 

We have implemented the above algorithms in Sage [24] for curves defined over 
Q p . In doing so, we made the following observations. 

4.1 Precision estimates 

For a tiny integral, the precision of the result depends on the truncation of 
the power series computed. Here is the analysis for a non- Weierstrass disc; the 
analysis for a Weierstrass disc, using a different local interpolation, is similar. 
(For points over ramified extensions, one must also account for the ramification 
index in the bound, but it should be clear from the proof how this is done.) 

Proposition 18. Let jp uj be a tiny integral in a non- Weierstrass residue disc, 
with P, Q defined over an unramified extension of K and accurate to n digits of 
precision. Let (x(t),y(t)) be the local interpolation between P and Q defined by 

x{t) = x(P)(l -t) + x(Q)t = x(P) + t(x(Q) - x(Pj) 

y(t) = y/fWfi). 



Let ui = g(x, y)dx be a differential of the second kind such that h(t) = g(x(t), y{t)) 
belongs to 0[\t]}. If we truncate h(t) modulo t m , then the computed value of the 
integral u> will be correct to min{n, m + 1 — [log p (r7j + l)j } digits of (absolute) 
precision. 

Proof. Let t' — t(x(Q) — x(P)). As P, Q are in the same residue disc and are 
defined over an unramified extension of K, we have v p (x(Q) — x(Pj) > 1. If we 
expand g(x(t'),y(t')) = X)to Ci (^')% then by hypothesis Ci e O. Thus 

rQ rQ 

J u = J g{x,y)dx 

= f g(x(t),y(t))dx(t) 
Jo 

x(Q)-x(P) 

g(x(t'),y(t'))dt' 

r x(Q)-x(P) oo 
= / Y^dit'fdt' 

OO 

= E-tt(-W)-^)) 1+1 - 

4=0 

The effect of omitting Ci(t') 1 from the expansion of g(x(t'),y(t')) for some i > m 
is to change the final sum by a quantity of valuation at least i + 1— \ \og p (i + l)\ > 
m+1— Llog p (m+l)J . The effect of the ambiguity in P and Q is that the computed 
value of (x(Q) — x{P)) l+1 differs from the true value by a quantity of valuation 
at least i + 1 — [\og p (i + 1)J + n — 1 > n. 

For Coleman integrals between different residue discs, which we may assume 
are non-Weierstrass thanks to § 3.4, one must first account for the precision loss 
in Algorithm 10. According to [17, Lemmas 2,3] and the erratum to [17] (or 
[15]), working to precision p N in Algorithm 10 produces the fi,Mij accurately 
modulo p N ~ n for n = 1 + [\og p max{ N, 2g + 1}J . 

We must then take into account the objects involved in the linear system (7), 
as follows. 

Proposition 19. Let fpUl be a Coleman integral, with lu a differential of the 
second kind and with P, Q in non- Weierstrass residue discs, defined over an 
unramified extension of Q p , and accurate to n digits of precision. Let Frob be 
the matrix of the action of Frobenius on the basis differentials. Set B = Frob* —I, 
and let m — v p (det(B)). Then the computed value of the integral fpU) will be 
accurate to n — max{m, [log p nj } digits of precision. 

Proof. By the linear system (7), the Coleman integral is expressed in terms 
of tiny integrals, integrals of exact forms evaluated at points, and a matrix 
inversion. Suppose that the entries of B = Frob* —I are computed to precision 




n. Then taking B~ x , we have to divide by det(B), which lowers the precision by 
m = v p (det(B)). By Proposition 18, computing tiny integrals (with the series 
expansions truncated modulo gives a result precise up to n— [\og p n\ digits. 

Thus the value of the integral J p lo will be correct ton- maxjm, [\og p n\ } digits 
of precision. 

4.2 Complexity analysis 

We assume that asymptotically fast integer and polynomial multiplication al- 
gorithms are used; specifically addition, subtraction, multiplication, and divi- 
sion take O(logiV) bit operations in Z/iVZ and 0(n) basering operations in 
R[x]/x n R[x]. In particular, this allows arithmetic operations in Q p to n (rela- 
tive) digits of precision, hereafter called field operations, in time 0(n logp). Using 
Newton iteration, both square roots and the Teichmuller character can be com- 
puted to n digits of precision using O(logn) arithmetic operations. (We again 
consider only points in non-Weierstrass discs defined over unramified fields.) 

Proposition 20. Let jp lo be a Coleman integral on a curve of genus g over Q p; 

with lo = df u + Y^!Li 1 c i w i a differential of the second kind and with P, Q in non- 
Weierstrass residue discs, defined overQ p , and accurate to n digits of precision. 
Let Frob be the matrix of the action of Frobenius on the basis differentials, and 
let m = v p (det(Frob* -/)). Let F(n) be the running time of evaluating f u at P 
and Q to n digits of precision. The value of the integral Jp lo can be computed 

to n — max{m, [log p n\ } digits of precision in time F(n) + 0(pn 2 g 2 + g 3 n\ogp). 
(Over a degree N unramified extension of Q p , the analysis is the same with the 
runtime multiplied by a factor of N.) 

Proof. An essential input to the algorithm is the matrix of the action of Frobe- 
nius, which can be computed by Kedlaya's algorithm to n digits of precision 
in running time 0(pn 2 g 2 ). Inverting the resulting matrix can be (naively) done 
with 0(g 3 ) arithmetic operations in Q p . It remains to be shown that no other 
step exceeds these running times. For the tiny integral on the first basis differen- 
tial, the power series x(t)/y(t) = x{t)f(x{t))~ 1 / 2 can be computed modulo t 11 ^ 1 
using Newton iteration, requiring O(nlogn) field operations. Each other basis 
differential can be computed from the first by multiplication by the linear poly- 
nomial x(t) and the definite integral evaluated with 0(n) field operations, for a 
total of 0(gn 2 ) bit operations. Computing <j)(P) and </>(Q) to n digits of preci- 
sion is cheap; directly using the formula in Algorithm 10 uses 0(g + logp) field 
operations. The last potentially significant step is computing and evaluating the 
/, at each P and/or Q. The coefficients of the fi can be read off in the reduction 
phase of Kedlaya's algorithm, and have 0(png) terms each. Evaluating (or even 
recording) all g of these forms takes 0(png 2 ) field operations, or 0(pn 2 g 2 ) bit 
operations, which is proportional to the cost of doing the reduction. 



4.3 Numerical examples 

Here are some sample computations made using our Sage implementation. Ad- 
ditional examples will appear in the first author's upcoming PhD thesis. 

Example 21. Leprevost [21] showed that the divisor (1, —1) — oo+ on the genus 
2 curve y 2 = (2x — l)(2x 5 — x 4 — Ax 2 + 8a; — 4) over Q is torsion of order 29. 
Consequently, the integrals of holomorphic differentials against this divisor must 
vanish. We may observe this vanishing numerically, as follows. Let 

~ 9 c 33 a 3 o 3 n 1 1 

C:y =x + w x + r + r ~r + T Q 

be the pullback of Leprevost's curve by the linear fractional transformation x i-> 
(1 — 2x)/(2x) taking oc to 1/2. The original points (1,-1), oo + correspond to 
the points P = (—1,1), Q = (0, j) on C. The curve C has good reduction at 
p = 11, and we compute 



ujq= I cji = 0(11 6 ), / uj 2 = 7-ll + 6-ll 2 + 3-ll 3 + ll 4 + 5-ll 5 + 0(ll 6 ), 



consistent with the fact that Q — P is torsion and cj ,^i are holomorphic but 
102 is not. 

Example 22. We give an example arising from the Chabauty method, taken from 
[23, § 8.1]. Let X be the curve 



whose Jacobian has Mordell-Weil rank 1. The curve X has good reduction at 7, 
and 



By [23, Theorem 5.3(2)], we know \X{Q)\ < 10. However, we can find 10 rational 
points on X: the six rational Weierstrass points, and the points (3, ±6), (10, ±120) 
Hence |X(Q)| = 10. 

Since the Chabauty condition holds, there must exist a holomorphic differ- 
ential to for which J® uj = for all Q e X(Q). We can find such a differential 
by taking Q to be one of the rational non- Weierstrass points, then computing 
a := luq, b := lj\ and setting u> = buiQ — auj\. For Q — (3, 6), we obtain 



Remark 23. It is worth pointing out some facts not exposed by Example 22. For 
instance, since to is already determined by a single rational non- Weierstrass point, 
we could have used it instead of a brute-force seach to find other rational points. 





y 2 = x(x - l)(x - 2)(x - 5)(x - 6), 



X(¥ 7 ) = {(0, 0), (1, 0), (2, 0), (5, 0), (6, 0), (3, 6), (3, -6), oo}. 




More seriously, in other examples, the integral w may vanish at a point defined 
over a number field which has a rational multiple in the Jacobian. Such points 
may be difficult to find by brute-force search; it may be easier to reconstruct 
them from p-adic approximations, obtained by writing f^ui as a function of a 
linear parameter of a residue disc, then finding the zeroes of that function. 

5 Future directions 

Here are some potential extensions of our computation of Coleman integrals. 
5.1 Iterated integrals 

Coleman's theory of integration is not limited to single integrals; it gives rise 
to an entire class of locally analytic functions, the Coleman functions, on which 
antidifferentiation is well-defined. In other words, one can define integrals 



These appear in several applications of Coleman integration, e.g., p-adic regula- 
tors in if-theory, and the nonabclian Chabauty method. 

As in the case of a single integral, one can use Frobenius equivariance to 
compute iterated Coleman integrals on hyperelliptic curves. One obtains a linear 
system expressing all n-fold integrals of basis differentials in terms of lower order 
integrals. Note that the number of such n-fold integrals is (2c/)™, so this is only 
feasible for small n. The cases n < 4 are already useful for applications, but 
ideas for reducing the combinatorial explosion for larger n would also be of 
interest. (One must be slightly careful in dealing with Weierstrass residue discs; 
see Remark 17.) 

We have made some limited experiments with double Coleman integrals in 
Sage. The Fubini identity 




turns out to be a useful consistency check for both single and double integrals. 
5.2 Beyond hyperelliptic curves 




which behave formally like iterated path integrals 




It should be possible to convert other algorithms for computing Frobenius ac- 
tions on dc Rham cohomology, for various classes of curves, into algorithms for 



computing Coleman integrals on such curves. Candidate algorithms include the 
adaptation of Kedlaya's algorithm to superelliptic curves by Gaudry and Giirel 
[14], or the general algorithm for nondegenerate curves due to Castryck, Denef, 
and Vercauteren [7] . It should also be possible to compute Coleman integrals us- 
ing Frobenius structures on Picard-Fuchs (Gauss-Manin) connections, extending 
Lauder's deformation method for computing Frobenius matrices [20]. 

5.3 Heights after Harvey 

We noted earlier that our algorithms for Coleman integration over Q p have linear 
runtime dependence on the prime p, arising from the corresponding dependence 
in Kedlaya's algorithm. In [15], Harvey gives a variant of Kedlaya's algorithm 
with only square-root dependence on p (but somewhat worse dependence on 
other parameters), by reorganizing the computation so that the dominant step 
is finding the p-th term of a linear matrix recurrence whose coefficients are 
polynomials in the sequence index. Harvey demonstrates the practicality of his 
algorithm for primes greater than 2 50 , which may have some relevance in cryp- 
tography for finding curves of low genus with nearly prime Jacobian orders. 

It should be possible to use similar ideas to obtain square-root dependence 
on p for Coleman integration, by constructing a recurrence that computes not 
just the entries of the Frobenius matrix but also the values fi(P) and fi(Q). 
However, this is presently a purely theoretical question, as we do not know of 
any applications of Coleman integration for very large p. 
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